In [1]:
pwd
Out[1]:
'C:\\Users\\alpha\\Documents\\TestingGeopandas'
In [1]:
cd C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents
C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents
In [2]:
import geopandas as gpd
import matplotlib.pyplot as plt
import pickle
import pandas as pd
import numpy as np
import calendar

import matplotlib.cm as cm
import matplotlib.colors as cls
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
In [3]:
#latex
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'custom'
matplotlib.rcParams['mathtext.rm'] = 'Times New Roman'
matplotlib.rcParams['mathtext.it'] = 'Times New Roman:italic'
matplotlib.rcParams['mathtext.bf'] = 'Times New Roman:bold'
In [17]:
from scipy.stats import boxcox, probplot
In [4]:
bios = gpd.read_file('biomes')
In [17]:
bios.columns
Out[17]:
Index(['OBJECTID', 'BIOME_NUM', 'SUM_2001_J', 'SUM_2001_F', 'SUM_2001_M',
       'SUM_2001_A', 'SUM_2001_1', 'SUM_2001_2', 'SUM_2001_3', 'SUM_2001_4',
       ...
       'SUM_2017_3', 'SUM_2017_4', 'SUM_2017_S', 'SUM_2017_O', 'SUM_2017_N',
       'SUM_2017_D', 'SUM_Area', 'Shape_Leng', 'Shape_Area', 'geometry'],
      dtype='object', length=210)
In [5]:
area = bios['SUM_Area']*1e-6
In [6]:
data_columns =bios.columns[2:2+17*12]
In [7]:
inter_cols = {}
for yr in range(17):
    inter_cols[yr] = data_columns[12*yr:12*yr +12]
In [8]:
intra_cols = {}
for mo in range(12):
    intra_cols[mo] = [data_columns[mo + yr *12  ] for yr in range(17)]
In [9]:
for yr in range(17):
    inter = bios[  inter_cols[yr]].sum(axis=1)/area
    label = 'inter_' + str(yr +2001)
    bios= bios.assign(label=inter)
    bios = bios.rename(columns={'label': label})
In [10]:
inter_columns = bios.columns[-17:]
In [11]:
inter_max = bios[inter_columns].max().max()
In [12]:
for mo in range(12):
    intra = bios[ intra_cols[mo]].sum(axis=1)/(17*area)
    label = 'intra_{:02}'.format(mo)
    bios = bios.assign(label=intra)
    bios = bios.rename(columns={'label': label})
In [13]:
intra_columns = bios.columns[-12:]
In [14]:
intra_max = bios[intra_columns].max().max()
In [36]:
bios.plot(column='intra_01', cmap='OrRd')
Out[36]:
<matplotlib.axes._subplots.AxesSubplot at 0x1e005d879e8>
In [15]:
def draw_map(df,column,title=None,vmin=0,vmax=1,name='name',dpi=100):
    fig, ax = plt.subplots(figsize=(12, 8))

    ax.set_aspect('equal')
    ax = df.plot(ax =ax, column=column, cmap='OrRd',vmin=vmin, vmax=vmax)
    ax.set_axis_off()
    
    norm = cls.Normalize(vmin, vmax)
    cmmapable = cm.ScalarMappable(norm, 'OrRd')
    
    cmmapable._A = []
    cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
    cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal") 

    cb.set_label(title, fontsize=20, family='Times New Roman')
    cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')

    plt.show()
    fig.savefig(name, dpi=dpi, bbox_inches='tight')

Inter-anual

In [38]:
%%time
for yr in range(17):
    column = 'inter_{}'.format(yr +2001)
    yr_name = str(yr +2001)
    title = r' Burned Area Pixels/$\mathrm{Km}^{2}$ by biome 2001-2017 ' +yr_name
    #title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
    #title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
    name = 'inter_biomes_{}'.format(yr +2001)
    draw_map(bios,column,title=title,vmin=0,vmax=inter_max,name=name)
    
Wall time: 5min 9s

Intra-anual

In [39]:
%%time
for mo in range(12):
    column = 'intra_{:02}'.format(mo)
    mo_name = calendar.month_name[mo + 1]
    title = r' Burned Area Pixels/$\mathrm{Km}^{2}$ by biome 2001-2017 ' + mo_name
    #title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
    #title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
    name = 'intra_biome_{:02}'.format(mo)
    draw_map(bios,column,title=title,vmin=0,vmax=intra_max,name=name)
Wall time: 3min 38s
In [42]:
bios[data_columns].sum(axis=1).sum()
Out[42]:
333598266.0
In [38]:
bios[data_columns].sum(axis=1)
Out[38]:
0           160.0
1      22274430.0
2       6778023.0
3        632711.0
4       6438595.0
5       1068629.0
6       4943428.0
7     238412159.0
8      11641973.0
9      14678367.0
10      1821057.0
11       465112.0
12      1632521.0
13     22519178.0
14       291923.0
dtype: float64
In [39]:
total = bios[data_columns].sum(axis=1)/17
In [40]:
total = total/area
In [41]:
bios= bios.assign(total=total)
In [46]:
column = 'total'
title = r' Burned Area Pixels/$\mathrm{Km}^{2}$ by biome 2001-2017 ' 
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'total_biomes'
draw_map(bios,column,title=title,vmin=0,vmax=total.max(),name=name)
In [47]:
column = 'total'
title = r' Burned Area Pixels/$Km^2$ by biome 2001-2017 ' 
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'total_biomes'
draw_map(bios,column,title=title,vmin=0,vmax=total.max(),name=name)
In [37]:
bios.total
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-37-09dac298e0da> in <module>()
----> 1 bios.total

C:\ProgramData\Anaconda3\envs\geopandas\lib\site-packages\pandas\core\generic.py in __getattr__(self, name)
   3612             if name in self._info_axis:
   3613                 return self[name]
-> 3614             return object.__getattribute__(self, name)
   3615 
   3616     def __setattr__(self, name, value):

AttributeError: 'GeoDataFrame' object has no attribute 'total'

Transformação Box-Cox

Inter-anual

In [18]:
lambdas = []
for yr in range(17):
    column= 'inter_{}'.format(2001 + yr)
    test_column, lbd = boxcox(bios[column] + 1e-25)
    lambdas.append(lbd)
lmbda=  sum(lambdas)/len(lambdas)
In [28]:
probplot(bios.inter_2017, plot=plt.gca());
In [20]:
probplot(test_column, plot=plt.gca());
In [23]:
lambdas
Out[23]:
[0.13985774858096542,
 0.20089681700197054,
 0.22312560927481467,
 0.1931761084130136,
 0.19599113319388284,
 0.20170717373301886,
 0.1959446551880218,
 0.20356773887009535,
 0.1742525330097433,
 0.1484901847211567,
 0.20684338419385906,
 0.2105950201665177,
 0.19555673704126555,
 0.14559862462925618,
 0.14543876619549959,
 0.19065394645742922,
 0.15070303733294083]
In [22]:
maxs = []
mins = []
for yr in range(17):
    column= 'inter_{}'.format(2001 + yr)
    test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
    maxs.append(test_column.max())
    mins.append(test_column.min())
In [24]:
box_min = min(mins)
box_delta =max(maxs) - box_min
In [25]:
for yr in range(17):
    column = 'inter_{}'.format(yr +2001)
    test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
    bios = bios.assign(test_column =test_column)
    yr_name = str(yr +2001)
    #"Burnt area pixels (#/km2) by biome in 2001
    title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome in ' +yr_name
    #title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
    #title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
    name = 'inter_biomes_index_{}'.format(yr +2001)
    draw_map(bios,'test_column',title=title,name=name)
C:\ProgramData\Anaconda3\envs\geopandas\lib\site-packages\matplotlib\font_manager.py:1339: UserWarning: findfont: Font family ['cursive'] not found. Falling back to DejaVu Sans
  (prop.get_family(), self.defaultFamily[fontext]))

Intra_anual

In [26]:
lambdas = []
for mo in range(12):
    column= 'intra_{:02}'.format(mo)
    test_column, lbd = boxcox(bios[column] + 1e-25)
    lambdas.append(lbd)
lmbda=  sum(lambdas)/len(lambdas)
In [27]:
lambdas
Out[27]:
[0.10605867201889993,
 0.09902942988287809,
 0.23832097659803791,
 0.15385833747596858,
 0.24977133293881681,
 0.2053018809757906,
 0.1845236015405936,
 0.19985762659223053,
 0.18808293768516446,
 0.12491423972870173,
 0.10850772664968505,
 0.09891453939720411]
In [29]:
probplot(bios.intra_11, plot=plt.gca());
In [30]:
probplot(test_column, plot=plt.gca());
In [31]:
maxs = []
mins = []
for mo in range(12):
    column= column= 'intra_{:02}'.format(mo)
    test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
    maxs.append(test_column.max())
    mins.append(test_column.min())
In [32]:
box_min = min(mins)
box_delta =max(maxs) - box_min
In [34]:
for mo in range(12):
    column = 'intra_{:02}'.format(mo)
    test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
    bios= bios.assign(test_column =test_column)
    mo_name = calendar.month_name[mo + 1]
    #"Burnt area pixels (#/km2) by biome in Janury (2001-2017)"
    title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome in ' +mo_name + ' (2001-2017)'
    #title = r' Burned Area Pixels/$\mathrm{Km}^{2}$ by climate 2001-2017 ' + mo_name
    #title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
    #title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
    name = 'intra_biomes_index_{:02}'.format(mo)
    draw_map(bios,'test_column',title=title,name=name)

Total

In [42]:
column = 'total'
test_column, lbd = boxcox(bios[column] + 1e-25)
test_column = (test_column - test_column.min())/(test_column.max() - test_column.min())
bios = bios.assign(test_column =test_column)
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome in (2001-2017)' 
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'total_biomes_index'
draw_map(bios,'test_column',title=title,name=name)
In [36]:
bios.columns
Out[36]:
Index(['OBJECTID', 'BIOME_NUM', 'SUM_2001_J', 'SUM_2001_F', 'SUM_2001_M',
       'SUM_2001_A', 'SUM_2001_1', 'SUM_2001_2', 'SUM_2001_3', 'SUM_2001_4',
       ...
       'intra_03', 'intra_04', 'intra_05', 'intra_06', 'intra_07', 'intra_08',
       'intra_09', 'intra_10', 'intra_11', 'test_column'],
      dtype='object', length=240)